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We study properties of the human electrocardiogram under 
the working hypothesis that fluctuations beyond the regular 
structure of single cardiac cycles are unpredictable. On this 
background we discuss the possibility to use the phase space 
embedding method for this kind of signal. In particular, the 
specific nature of the stochastic or high dimensional compo- 
nent allows to use phase space embeddings for certain signal 
processing tasks. As practical applications, we discuss noise 
filtering, fetal ECG extraction, and the automatic detection 
of clinically relevant features. The main purpose of the pa- 
per is to connect results of embedding theory which had not 
been previously applied in practise, and practical applications 
which had not yet been justified theoretically. 



I. INTRODUCTION 

In nonlinear time series analysis, methods developed 
in nonlinear dynamics are applied to time series data in 
order to capture as much of the underlying structure as 
possible. One central question hereby is whether the dy- 
namics of the system can be reconstructed on the basis of 
the data given. Given the case that the considered system 
is deterministic, the dynamics unambiguously evolves one 
state in phase space into another. Therefore, the recon- 
struction of phase space is a fundamental problem which 
plays the key role in many applications: as long as the 
state space reconstruction remains unjustified all of the 
consequent analysis may be wrong. From this point of 
view, we discuss the phase space embedding of realistic 
signals in the example of the normal human electrocar- 
diogram. An introduction to nonlinear dynamics is found 
in |l],^| . An account of nonlinear time series methods is 
given by [g. 

The electrocardiogram (ECG) is one of the most 
prominent clinical tools to monitor the activity of the 
heart. In order to record an electrocardiographical sig- 
nal, metal electrodes are placed on the patient's chest 
wall and extremities (for details see e.g. Q). The poten- 
tials are generated by the atrial and ventricular muscle 
fibers. Due to the placement of the electrodes on the 
skin at some distance from the heart, the signals mea- 
sured correspond to action potentials which are averaged 
over large regions of tissue. The spreading of the elec- 
trical activity over the cardiac muscle is controlled by 
the conduction system of the heart. Atrial and ventric- 
ular contraction and relaxation, respectively, correspond 
to characteristic ECG-waves which are traditionally la- 
beled in alphabetic order beginning with the letter 'P'. 
As long as the average heart rate does not change dra- 
matically, a non-pathological ECG shows a nearly peri- 




FIG. 1. Two different cycles of the same ECG have been 
aligned (indicated by the arrow) right before the beginning of 
the P-waves. The stochastic components which lead to devi- 
ations from a pure limit cycle behavior manifest themselves 
in different cycle lengths and a slight variation of shape. 



odic structure which is due to the continuous generation 
of action potentials and the fixed pattern according to 
which the electrical activity spreads out over the cardiac 
muscle. However, apart from the deterministic structure 
some kind of variability can be found. On the one hand, 
the length of the time interval between successive beats 
fluctuates - to a certain extend - around the mean heart 
rate. (In the medical literature the inter beat-interval 
is often called RR-interval which is defined as the time 
between two consecutive R-waves). On the other hand, 
a slight variation of the cycle shape and amplitude can 
be observed. Both kinds of fluctuations can be gath- 
ered from Fig. [l] where two different cycles of the same 
ECG sequence are shown. As one can see in Figure ||, 
part of the variation of the heart rate (middle trace) of 
a human at rest can be connected with the breath cycle 
(upper trace). (Data from the Santa Fe Institutte time 
series contest || , see Ref . Q ) . The lower trace contains a 
randomized sequence that has the same autocorrelation 
structure, and the same cross correlation to the breath 
rate as the middle trace. For details as the generation of 
these surrogate data, and further references see Ref. [Q. 
This random surrogate data explain most but not all of 
the underlying structure by linear correlations. One ex- 
planation for the remaining structure might be the pres- 
ence of a high dimensional or nonlinear stochastic com- 
ponent. In this paper we adopt as a working hypothesis 
that the fluctuations in the instantaneous heart rate are 
effectively unpredictable. 

The Fourier domain is inappropriate to capture the 
structure described above, since the variation in cycle 
length leads to a dominating broad band component. 
An alternative would be to formulate the dynamics of 
the ECG cycle by a stochastically driven model in a low 
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FIG. 2. Simultaneous measurements of breath and heart 
rates, upper and middle trace. The lower trace contains a 
randomized sequence preserving the autocorrelation structure 
and the cross correlation to the fixed breath rate series. The 
surrogate data mimic much of the structure contained in the 
data, but not all. 



dimensional phase space. Since this space is not fully 
accessible by measurements, we propose to use the delay 
reconstruction technique ||] as a convenient tool to reveal 
both the regular and stochastic aspects of the electrocar- 
diogram. Hereby, one attempts to reconstruct the state 
variables £ of the system represented in phase space on 
the basis of the single lead measurement. Phase space 
reconstruction has been developed for purely determin- 
istic systems. We do not want to make this assumption 
here. We will argue, however, why the technique can be 
appropriate in this case despite the stochastic component 
of the signal. 



A. Delay reconstruction of electrocardiograms 

In many cases of practical interest - like in the ECG 
- it is not possible to measure the state variables of a 
system directly. Instead, the measuring procedure yields 
some value x = when the system is in state 

Here, if is a measurement function which in general de- 
pends on the state variables in a nonlinear way. The 
time evolution of the state of the system results in a 
scalar time series Xi,X2,X3,... . In order to recon- 
struct the underlying dynamics in phase space, delay em- 
bedding techniques are commonly used. Delay vectors 
x„ = (x n ,x n -i,x n _2i, ■ ■ ■ ,£„_(d_i);) are a convenient 
method to transform the scalar time series into d dimen- 
sional vectors. Herein, d corresponds to the embedding 
dimension while / is the lag between the time series ele- 
ments. In the case of purely deterministic systems, the 
embedding theorem by Takens |8| and its generalization 
by Sauer et al. |J assure the equivalence between the 
reconstruction space and the original phase space of the 
system under fairly broad conditions. Figure shows a 
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FIG. 3. Two dimensional delay reconstruction with delay 
time 12 ms using 5000 data points of the whole ECG sequence 
from which two cycles are shown in Fig. [j]. 



two dimensional reconstruction of the first 5000 time se- 
ries elements (20s) of the ECG sequence underlying Fig. |l| 
by using delay coordinates at a lag of 12 ms. 

Whilst for the signal processing applications discussed 
later usually much higher embedding dimensions are 
needed, even the two dimensional representation illus- 
trates the underlying phase space structure succinctly. 
It is remarkable that the signal structure in phase space 
is close to a manifold (which will be exploited for noise 
reduction purposes discussed later). This suggests that 
the underlying generating process is basically determin- 
istic which manifests itself in regular cycles, while the 
stochastic degrees of freedom lead to deviations from the 
basic cycle. While the fluctuations in cycle shape are in- 
dicated in the delay reconstruction, the stochastic varia- 
tion of the time between cycles cannot be seen directly. 
This can be easily understood in the context of the fol- 
lowing simplifying model. Let us suppose that the onset 
of a new cycle is started by an external trigger signal 
and that the RR-interval time fluctuates stochastically 
around the mean heart rate. Let us further assume that 
the ECG signal rests on the base line after a cycle is 
finished, until a new trigger signal starts the next depo- 
larisation. Since the constant base line corresponds to a 
single point in phase space, the stochastic variation of the 
cycle length cannot be revealed using delay vectors, pro- 
vided they span a time shorter than the interbeat times. 
In order to reconstruct the underlying dynamics fully, it 
is necessary to use some additional information - like the 
trigger signal - describing the onset of the next cycle. 

To conclude this introduction, the delay reconstruction 
technique can help to reveal aspects of the nature of the 
underlying process, but it is not clear a priori that the 
dynamics of the system can be reconstructed by using 
delay coordinates in the case of stochastic forcing, even 
if the embedding dimension is chosen sufficiently large. 
Fortunately, the embedding theorems proven recently by 
Stark and coworkers (in the following we refer to 

these theorems by the initials of the authors: SBDO) 
generalize the well known Takens embedding theorem B 
and provide results for the phase space reconstruction of 
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stochastically forced dynamical systems. On the basis 
of the SBDO-theorem it is possible to reconstruct the 
higher dimensional phase space from the scalar time se- 
ries data by using delay coordinates when the sequence of 
stochastic influences that act on the system is known. Of 
course, in the case of the ECG, we don't know the driv- 
ing process. However, we will argue that the necessary 
information can be recovered from the signal a posteriori. 



II. EMBEDDING OF STOCHASTICALLY 
DRIVEN SIGNALS 

In this section we give an overview over the funda- 
mental notions introduced in the papers by Stark and 
coworkers [[lO 11 in order to state the key points of the 
essential embedding theorem for stochastic systems. We 
try to simplify the notation used as far as possible. The 
statement of the theorem as well as the formalism in- 
troduced are illustrated in the examples discussed subse- 
quently. 

In order to describe stochastic forcing, Stark considers 
skew product systems (to avoid confusion we denote the 
state variables of the system by £ while x labels the delay 
vectors build from the scalar time series measurements in 
order to reconstruct the phase space of the system) : 



€„+i = f(£„,n») 

The dynamical system f is stochastically driven by the 
bi-infinite sequence = . . . lu^2^-i<-^q'-^i^ i 2 ■ ■ ■, where Q 
is build up by symbols u>i that influence f . The time evo- 
lution of the stochastic sequence f2 is realized by the map 
a shifting the elements of f2 to the left by one position. In 
each time step, the dynamics f is affected only by the cen- 
tral element u>q of Q. It is natural to express the depen- 
dence of f on oj n by a sequence of maps f„(£) = f(£, w n ), 
so that the dynamical system becomes: 



(1) 



In this framework it is possible to modify the underly- 
ing dynamics in each time step. This formalism will be 
illuminated in an example given below where the phase 
space reconstruction of a circular motion with fluctuating 
radius is discussed. 

The central theorem of jll| faces the question whether 
it is possible to reconstruct the phase space of stochasti- 
cally driven systems only on the basis of the time se- 
ries. In analogy to the usual embedding techniques, 
one can define the d dimensional delay embedding map 
$ : M — > lZ d for skew product systems by 

MO = ( <P (O > V (Ol) ,...,¥> (Zn+d-l)) (2) 

(the usage of "prelay" rather than delay vectors is for 
technical reasons). Here, the index u> has to be kept 



since the future of £„ depends explicitly on the stochastic 
sequence f2 under consideration. This can be stressed by 
writing 

*«(0 = (v(0>¥>(k(0),"- 

...,^(f n+d _ 2 (...f n+1 (f„(0) •••))) (3) 

where we have made use of eqn. ([!]). Essentially, the 
main Theorem (No 3.5 JllJ] and No 7 fio) respectively) 
says that if d > 2do + 1, then there is a residual set of 
dynamical systems f and measurement functions ip, so 
that the delay embedding map $ w yields an embedding 
for almost every stochastic sequence u>. For the detailed 
genericity conditions on f , ip and lo we refer to pp| , pj| . 

The crucial point here - compared with the embedding 
theorems by Takens (^] and Sauer et al. ^ — is that it is 
not possible to reconstruct the underlying dynamics fully 
on the basis of only the measured time series, even if the 
embedding dimension d is sufficiently large and the dy- 
namical system f as well as the measurement function ip 
are generic. In order to take the stochastic nature of the 
process into account, it is necessary to use additional in- 
formation describing the stochastic influences. To be pre- 
cise, the knowledge of the underlying stochastic sequence 
is essential for the definition of the delay embedding map 



As an illustrative example consider the map 
x n+ i \ r n+ i. 



Vn+l 



-Rr, 



(4) 
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where R^ is a rotation by 



Rrf. 
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and the radius r n is chosen at random between r m i n and 
^max while the rotation angle <f> is fixed and generically 
incommensurate with 2n. 

For simplicity, let the measurement function tp 
sure" the first coordinate: 



mea- 



X n 



The central question is now whether the scalar time se- 
ries xq, . . . ,x n , x n +i, ■ ■ ■ yields an appropriate basis to 
reconstruct the dynamics of the system. Fig. |] contains 
a schematic representation of the phase space structure 
underlying equation (^). Since the fluctuations of the ra- 
dius are bounded by r min and r max , the points generated 
according to Eq. (Q) are sprinkled over a disk. This disk 
corresponds to the compact manifold M in the theorems 
by Stark et al. 

For a moment consider the case that the sequence of 
radii would be fixed to one: r n = 1. In this specific 
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FIG. 4. Illustration of the phase space structure underlying 
equation (^). After many iterations, the dots evenly fill the 
disk area between r m i n and r max . 



case it is obvious that the two dimensional delay recon- 
struction map $(x„_i) = (x n ^i,x n ) - which is based 
on the measurement function specified above - yields an 
embedding for almost every value of (/). This is simply 
due to the fact that a uniform motion on a circle can 
easily be reconstructed by its hrst Cartesian coordinate 
since both coordinates are related by a constant phase 
shift. As far as a deterministic variation of the radius is 
involved, - according to Takens' theorem - at least three 
dimensional delay embedding maps are needed to avoid 
self intersections in reconstruction space. It can be seen 
quickly that a pure delay embedding technique in which 
the dimension takes into account only the degrees of free- 
dom of f is in general not sufficient to reconstruct the 
underlying dynamics in the case of stochastic forcing. In 
contrast to deterministic systems, due to the stochastic 
influences crossing points can appear between the tra- 
jectories on the manifold M as indicated in Fig. [|. In 
order to determine the continuation of the trajectories, 
it is necessary to use some information beyond the delay 
coordinates. Fortunately, as Stark points out in the theo- 
rem cited above, the additional knowledge of the driving 
stochastic sequence allows the full reconstruction of the 
underlying dynamics. Referring to map (Q), the knowl- 
edge of the sequence of radii r„ allows the continuation 
of the trajectories beyond the crossing points. 

Apart from the pure illustrating function regarding the 
basic idea underlying Stark's theorem, the model dis- 
cussed can help to motivate the embedding of more real- 
istic signals like ECG sequences. As indicated in fig. [|, 
the delay reconstruction of a normal human electrocar- 
diogram bears resemblance to a limit cycle structure with 
fluctuations. As pointed out in the introduction, the fluc- 
tuations affect both the cycle length and its shape. One 
natural way to describe the stochastic influences would be 
by trying to observe the driving input signal. However, 
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FIG. 5. Trajectories of a kicked, damped harmonic oscil- 
lator. Upper: true phase space. Lower: delay embedding. 
Kicks occur close in time. 



the variable parameters that influence the ECG wave- 
form are not usually accessible to measurements and in 
fact, are not fully known. This is a general problem with 
the approach behind the work of Stark et al., as well as 
that of Casdagli jl2|: regarding a system as an input- 
output device is only useful for the analysis of time series 
if both, the output and the input sequences are observed. 

One of the main ideas we want to put forth in this 
paper is that there are situations in which the input se- 
quence can to some extent be inferred from the observed 
data. In such a case, the embedding procedure may be 
applied successfully, on the theoretical foundation of the 
theorem discussed above. In the case of the ECG, the 
stochastic driving mainly affects the duration of the cy- 
cles, which is accessible a posteriori by measuring the RR 
intervals. If there is a variation of the ECG waveform it- 
self which can be captured by a few parameters, these 
can also be measured once the ECG has been recorded. 
The fact we use in doing this is that ECG recordings 
contain considerable redundancy. With regard to ECG 
signals this implies that the explicit specification of the 
underlying sequence Q is not strictly necessary to make 
use of the embedding technique. 
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FIG. 6. Same as Fig. ^, but kicks are well separated in time 
such that the system can relax between kicks. 



Let us study one more toy model to illustrate how 
the lacking information about the input sequence may 
be contained in a signal from an input-output system. 
Consider a damped harmonic oscillator 



x + x + x = a(t) , 



(5) 



where the driving term is zero except for kicks of ran- 
dom strength at times U such that the inter-beat inter- 
vals, U — ti-i are random in the interval [p, q\. Figures || 
and |^ show two trajectories of such a system with dif- 
ferent choices of the inter-beat time interval [p, q] . On 
discretisation, the kicks are realized by finite jumps by a 
random amount in the interval [0,1]. The upper panel 
in each of these figures shows the true phase space, while 
in the lower panels a delay representation is used with a 
delay of one time unit. In Fig. ^, beats are initiated 
with a time separation of p = to q — T/2 where 
T = 47r/\/3 ~ 7.26 time units is the period of oscilla- 
tion. This does not allow the system to relax sufficiently 
between kicks in order to form characteristic structure in 
phase space. Consequently, an embedding provides no 
clear picture. In Fig. fl, no kicks were closer in time than 
p = T, the maximal separation being q = 3T. The inter- 
beat parts of the trajectory are distinct because of the 



different kick strength, but since this is the only fluctu- 
ating parameter, they are essentially restricted to a two- 
dimensional manifold. This manifold is preserved under 
time delay embedding although neither the sequence of 
beat times nor the beat amplitudes are used explicitly. 
The randomness prevails in the indeterminacy at the ori- 
gin as to when the next beat will occur. 



III. PREDICTING AND PROJECTING ECG 
SIGNALS 

While the structure of a normal human ECG is well 
predictable on a scale considerably shorter than one cy- 
cle, fluctuations in cycle length limit the accuracy at the 
onset of the QRS-complexes. Therefore, large deviations 
can be found between the original ECG data and the 
one step predictions generated as soon as a new QRS- 
complex appears, whilst the difference between both sig- 
nals is usually comparable to the magnitude of the noise 
level. Fig. ^ illustrates these deviations with the ECG 
sequence plotted. The ansatz chosen makes explicit use 
of the phase space reconstruction by delay vectors x„. 
The ECG data considered are highly oversampled so that 
x n +\ is close to x n up to some corrections which are mod- 
eled by radial basis functions: 



K 

c n+ i = x n +^2a k exp (^-b [x„ - c k \ 



C 



(6) 



k=l 



F(x„) 



The constant C allows to handle a zero offset directly 
instead by fitting via radial basis functions, which would 
affect the quality of the parameters a k , b and Cfc. In order 
to place the centers c k of the basis functions, a grid-based 
algorithm was used. The result plotted in Fig. ^ was 
achieved by aspiring a uniform density of centers over the 
phase space regions occupied by the reconstruction vec- 
tors and a total number of 10 basis functions. Increasing 
the number of basis functions leads to a better approxi- 
mation of the ECG signal and therefore a lower in-sample 
prediction error. In Fig. ^ the prediction error during the 
interbeat times is comparable to the magnitude of the 
original ECG data noise level. It is remarkable that the 
error of the one-step predictions is much larger at the on- 
set of the QRS-complexes so that the overall distribution 
is nonuniform. Therefore one might conclude that the 
onset of the next QRS-complex remains unpredictable. 

To reduce the uncertainty underlying the onset of the 
next QRS-complex one could think of using the length 
of the following RR-interval as an additional input sig- 
nal. It might be possible to lower the large prediction 
error at the beginning R-wave by this technique, since 
the moment of the spike-like deflection is now completely 
determined. The use of an additional input sequence 
which characterizes the stochastic degrees of freedom cor- 
responds to the fixing of the stochastic sequence f2 in the 
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FIG. 7. Deviations between the ECG sequence |13| (upger 
trace) and one step predictions generated according to Eq. 
The signals are measured in uncalibrated AD-units. 



formalism introduced by Stark. Since the prediction of 
ECG signals is not in itself of much practical relevance, 
we do not describe further details. The concrete imple- 
mentation is technically involved and would go far be- 
yond the scope of this paper. 

State space reconstruction techniques form an essen- 
tial tool for noise reduction purposes. In order to sepa- 
rate the signal from its noise components, basically two 
different techniques can be pursued. (I) On the one hand 
a global fit to the dynamical equations can be used for 
the reconstruction of a noiseless trajectory of the system, 
as it has been done in [Q, for example. In this context, 
state space reconstruction techniques are involved to cap- 
ture the underlying dynamics. For example in Eq. (j^) 
the reconstruction technique enters via the delay vec- 
tors x„. As an iterative method for reaching a trajec- 
tory which obeys the reconstructed dynamics, gradient 
descent methods have been investigated in |L4j. (II) On 
the other hand it is possible to exploit the local structure 
in the reconstruction space directly for noise reduction 
purposes without referring to the global underlying dy- 
namics. Methods based on this local projective approach 
have been proposed in jl5| , see Jl6| for a review. The idea 
is that if it is found empirically that the data points in 
reconstruction space are located close to a manifold, the 
error by projecting onto that manifold may be smaller 
than the error due to the noise. If the trajectory of the 
system lies on a low dimensional attractor, the projec- 
tion technique can indeed be shown to reduce noise. But 
also for non-deterministic signals, nonlinear noise reduc- 
tion techniques have been successfully applied for signal 
processing tasks like ECG noise reduction Jlj],[l8| or fe- 
tal ECG extraction |^9|,^o|. The basic reason why this 
works is the mechanism illustrated in fig. ^|. These sig- 
nals, in spite of being non-deterministic, nevertheless are 
near a low-dimensional manifold. Let us give the extrac- 
tion of the fetal ECG as an application of the delay em- 
bedding of ECG signals (fetal electrocardiography is the 
only method to monitor the cardiac activity of the fetus 
in a non invasive way) . Figure |8| contains a two dimen- 
sional delay representation of an abdominal ECG record- 
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FIG. 8. Delay reconstruction of the abdominal ECG 
recording (left panel) and the maternal component (right 
panel) . 
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FIG. 9. Fetal ECG extraction by locally linear projections 
in phase space. The abdominal recording (upper trace) con- 
tains a dominant maternal component and the fetal signal 
which is covered by noise. The first noise reduction step yields 
the reconstructed maternal signal (middle) which allows to 
separate the noisy fetal ECG from the original recording. A 
second projection can be applied to clean up the fetal com- 
ponent (lower trace). 



ing (left panel) and the reconstructed maternal manifold 
structure (right panel) by the projective noise reduction 
approach. Segments of the corresponding time series are 
shown in Fig. || (data by courtesy of J. F. Hofmeister 
pl[). The upper trace contains the abdominal signal 
(electrodes have been placed on the maternal abdomen 
to record most of the fetal heart activity) while the mid- 
dle trace shows the abdominal projection of the mater- 
nal electrocardiogram which has been extracted by noise 
reduction. The difference between the upper and mid- 
dle trace yields the noisy fetal component, which can be 
cleaned up in a second noise reduction step (lower trace) . 
The extraction of the fetal electrocardiogram is broadly 
described in |^2| , for a powerfull modification of the pro- 
jective algorithm which allows the real time extraction 
of the fetal electrocardiogram on a Laptop PC (Pentium 
processor at 133 MHz, Linux operating system, 250 Hz 
sampling rate) see pl| . 

As another practical application, phase space recon- 
struction techniques can be employed for the localization 
of fundamental ECG waves like P-, R- or T-waves. In- 
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stead of comparing the structures directly on the basis of 
the time series, similar ECG waves can be identified as 
adjacent curves in reconstruction space. For example the 
smaller loop to the lower left side of the origin in fig. || 
corresponds to the P-waves of the full ECG sequence par- 
tially plotted in fig. [I], while the larger framing loop is 
due to the QRS-complex. For the automatic detection of 
ECG cycle substructures, delay vectors are defined to be 
neighbours if their distance to the vector under consider- 
ation is smaller than a given size (for neighbour searching 
methods cf. e.g. Nearby trajectories can be identified 
in phase space by searching for neighbours of the recon- 
struction vectors of a specified structure. On the basis 
of the original time series, these trajectories correspond 
to segments which are similar to the specified one. By 
tagging a P-, R- or T-wave of a given electrocardiogram, 
previous or following waves of this type can be located 
automatically. 

To conclude, embedding techniques form a fundamen- 
tal starting point for many applications. Therefore, it is 
desirable to motivate the reconstructability of the phase 
space in the case of typical time series. We illuminated 
this issue by considering the example of the human elec- 
trocardiogram. Under the working hypothesis that fluc- 
tuations beyond the regular structure of the cardiac cycle 
are unpredictable we found the embedding theorems of 
Stark and coworkers [ ^0|JTl] ] to be useful to motivate the 
embedding of ECG data. Finally, we discussed useful 
applications which are based on phase space reconstruc- 
tions of the electrocardiogram in order to give examples 
of how embedding techniques can be employed for prac- 
tical tasks. 
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